Linear response theory in the continuum for deformed nuclei: 
Green's function vs. time-dependent Hartree-Fock with the absorbing-boundary 

condition 
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The continuum random-phase approximation is extended to the one applicable to deformed nuclei. 
We propose two different approaches. One is based on the use of the three dimensional (3D) Green's 
function and the other is the small-amplitude TDHF with the absorbing-boundary condition. Both 
methods are based on the 3D Cartesian grid representation and applicable to systems without any 
symmetry on nuclear shape. The accuracy and identity of these two methods are examined with the 
BKN interaction. Using the full Skyrme energy functional in the small-amplitude TDHF approach, 
we study the isovector giant dipole states in the continuum for 16 O and for even-even Be isotopes. 

PACS numbers: 21.60.Jz, 21.10.Pc, 27.20.+n 



I. INTRODUCTION 

Mean-field theories with effective interactions 0, S 
4] have been extensively used for systematic description 
of nuclear ground-state properties from light to heavy 
nuclei, including infinite nuclear matter. Nuclear mass, 
radius, density distribution, and deformation are the pri- 
mary target of the static effective mean-field theory 0. 
The concept of the nuclear mean-field theory is rather 
different from the Hartree-Fock theory in electronic sys- 
tems but is more close to the density functional theory. 
Especially, the Hartree-Fock (HF) with the zero-range 
Skyrme interaction results in an energy functional of lo- 
cal densities. A similar form of functional was obtained 
from the density-matrix expansion of energy functionals 
calculated with the microscopic nucleon-nucleon forces 
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Although the static mean-field calculations well repro- 
duce the bulk nuclear properties throughout the nuclear 
chart, it is necessary to go beyond the mean field to 
describe excited states and correlations associated with 
many kinds of collective motions. The generator coordi- 
nate method (GCM) dH is one of the standard methods 
to take account of the configuration mixing. The GCM 
based on the mean-field theory provides a unified descrip- 
tion of single-particle and collective nuclear dynamics. In 
practice, collective variables, q, are chosen from physi- 
cal intuition and are restricted to one dimension in most 
cases. For instance, in order to describe quadrupole exci- 
tations, the most common choice is the mass quadrupole 
moment, q = A(§(q)\r 2 Y2o\${q)} , where the single-Slater 
states | $(q)) are determined by the constrained Hartree- 
Fock(-Bogoliubov) calculation. This is a drawback of the 
GCM that one has to prepare, a priori, a set of states 

The time-dependent Hartree-Fock (TDHF) theory is a 
complementary method to the GCM. The system deter- 
mines its collective path for itself and the TDHF takes 
care of both collective and single-particle excitations. 



The TDHF is also known to produce the proper iner- 
tial parameters |1C| . because it is a dynamical theory 
to incorporate time-odd components in the wave func- 
tion. A drawback is its semiclassical nature. Namely, 
in order to calculate quantal quantities, such as eigenen- 
ergy and transition probability, one has to requantize ob- 
tained TDHF dynamics. Although it is a difficult task 
to requantize the TDHF orbitals in general > the per- 
turbative regime can be easily handled. The linear ap- 
proximation leads to the random-phase approximation 
(RPA) for the effective density-dependent forces, which 
is analogous to the time-dependent local-density approxi- 
mation in electronic systems |l2l 1 1 3| . Another advantage 
of TDHF is its ability of describing spreading width of 
collective motion induced by the interaction between par- 
ticles and time-dependent mean-field potential (one-body 
dissipation). The escape width can be also described by 
the TDHF but requires proper treatment of the contin- 
uum. In this paper, we propose a feasible method to treat 
the continuum in the real-space TDHF calculation. That 
is the absorbing-boundary condition (ABC) approach. 
We have already studied photoabsorption in molecules 
fl4| and nuclear breakup reaction fl5lll6j| with the simi- 
lar technique. Our earlier attempts for nuclear res pon se 
calculation have been reported in Refs. [rH ITsl Il9l |2C( . 

The Green's function method in the linear response 
proposed by Shlomo and Bertsch plj is a common way 
to treat the continuum boundary condition. It is usu- 
ally called "continuum RPA" in nuclear physics. The 
same idea was proposed later in the time-dependent den- 
sity functional theory (TDDFT) for calculations of pho- 
toresponse in rare-gas atoms [22 ] . The method has been 
widely applied to spherical (magic or semi-magic) nu- 
clei [SlillSIMlllllllialHaEllll, however, its 
application to deformed systems has not been done so 
far, because the explicit construction of Green's function 
is extremely difficult for deformed potential. We have 
recently proposed an iterative method to construct re- 
sponse functions for deformed systems with the proper 
boundary condition in the three-dimensional (3D) coor- 
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dinate space and studied molecular photoabsorption us- 
ing the TDDFT [HIIJ. There, the dynamical screen- 
ing effect in the continuum for a multi-center problem 
was a key issue for understanding the photoabsorption 
cross section at photon energies higher than the ioniza- 
tion potential. The same method is applicable to nuclear 
mean-field models that do not contain non-local densi- 
ties. In this respect, applications to the Skyrme energy 
functional is parallel to the TDDFT. In this paper, we 
extend a method of the continuum RPA to the one in 
the 3D coordinate space and apply it to deformed nuclei. 
We call this "3D continuum RPA" in this paper. This 
provides the exact treatment of the nucleonic continuum 
for deformed nuclei Q ■ The results can be used to check 
validity of the ABC approach. 

Another issue addressed in this paper is the self- 
consistent treatment in the continuum response calcu- 
lation. The nuclear energy functional is far more com- 
plicated than that of electronic TDDFT. The Skyrme 
functional is one of the simplest, since its non-local part 
is expressed by derivatives of local densities. Even so, the 
continuum RPA calculations so far neglect the spin-orbit 
and Coulomb residual particle-hole interactions, which 
violates the self-consistency with the HF field p9l . [35| . 
In addition, a time-reversal-odd (time-odd) part of den- 
sities, such as spin densities, are often omitted. Since 
the spin-orbit term in the time-even mean field is re- 
lated to spin-current terms in the time-odd mean field 
by the local gauge (Galilean) invariance |M H3, the 
neglect of spin density violates this symmetry. As far 
as we know, at present, there is no fully self-consistent 
Skyrme-HF-based continuum RPA calculations, even for 
spherical nuclei. We perform the small-amplitude TDHF 
calculation with the ABC (TDHF+ABC) in fully self- 
consistent manner for the giant dipole resonance in 16 O, 
and examine effect of residual interactions which have 
been neglected so far. In the time-dependent relativistic 
mean-field approach without the continuum, the small- 
amplitude real-time calculation has been attempted for 
spherical nuclei |38ll39l|. However, only a very short time 
period (3^4 HMeV" 1 ) was achieved, which prevents 
them from carrying out a quantitative analysis. See also 
recent papers [33, Ell and references therein for the 
present status of the self-consistent HF(B)+(Q)RPA cal- 
culations for spherical nuclei. It should be noted that, 
without the continuum boundary condition, there exist 
a few works of fully self-consistent RPA for deformed nu- 
clei calculated in the 3D coordinate space with the full 
Skyrme interaction [43LE4| . 

In recent developments of radioactive-ion-beam facili- 
ties, the Coulomb excitation and the inelastic scattering 
are becoming standard methods to investigate excited 
states in unstable nuclei. For weakly bound systems, the 
treatment of the continuum should be extremely impor- 
tant. Moreover, we know most of open-shell nuclei arc 
deformed. Collective modes of excitation in the particle 
continuum in deformed nuclei become the main interest 
in those studies. This has not been examined in a self- 



consistent manner so far. The present paper will provide 
us with general methods of linear response in the con- 
tinuum for systems whose energy functional is given by 
local one-body densities. 

The paper is organized as follows: Section ITT1 presents 
a method of extending the continuum RPA to deformed 
nuclei. In Sec. lIIII we present a real-time TDHF method 
using the absorbing boundary condition. Some illustra- 
tive examples show effect of the continuum and compari- 
son between these two methods in Sec. IIVI In Sec. 
we present numerical results of the small-amplitude 
TDHF+ABC calculation in real time using the Skyrme 
energy functional for giant dipole resonances. Effects of 
time-odd densities, El strengths in 16 0, and those in 
neutron-rich Be isotopes are discussed. The conclusion 
is summarized in Sec. I VII 

II. 3D CONTINUUM RPA 

For spherical systems, the continuum RPA is formu- 
lated in terms of the radial Green's function using a 
multipole expansion |2l|. Hereafter, we refer to this as 
"ID continuum RPA". In Ref. [l4|, we have presented a 
method to construct a Green's function in the 3D grid 
representation for a system without any spatial symme- 
try. In this section, we recapitulate the method of con- 
structing the response function. Spin and isospin indices 
are suppressed for simplicity and h — 1 is used. 

The HF Hamiltonian, h[p], is a functional of one-body 
density matrix 45| . In case of zero-range effective inter- 
actions, it is a functional of local one-body density p(r). 
The stationary condition is 

[M = o, (i) 

which defines the HF ground state density p = pq. Then, 
the TDHF equation with an external perturbation, 

ij t p(t) = [h[p] + V ext (t),p(t)}, (2) 

is linearized with respect to the density fluctuation, 

p(r,t) = po(r) + Sp{r,t). (3) 

This leads to the well-known RPA equation. The transi- 
tion density, 8p(r\ui), which is the Fourier transform of 
Sp(r, t), can be expressed as [22ll35j 

Spfow) = J J drU(r,v';u;)V e ^(r';co), (4) 

- J dr'n (r,r'; W ) 

x (v cxt (r';cj) + J dr"v(r',r")8p(r";tj))(,5) 

where n and n are the RPA and the independent- 
particle response function, respectively. The u(r, r') is 
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a residual interaction which is defined by 

S 2 E[p] 



v(r,r') 



Sp(r)6p(r') 



(6) 



Here, we assume that the total energy functional, E[p], 
is a function of local density p(r) only. The HF mean 
field is also local in the coordinate space. Assuming that 
a one-particle moment F(r) depends only on the spatial 
coordinates, the transition strength is obtained from the 
transition density, 



dB{uj;F) 
du> 



\(n\F\0)\ 2 6(c -E n ), (7) 

n 

--Im J J rfrdr'F(r)n(r,r';w) J F 1 (r / )(8) 
-- Im / drF(r)6p{r;u) . (9) 

J V^ = F 



In case of deformed nuclei, |0) and \n) are not eigen- 
states of total angular momentum operator. Thus, 
dB(u>; F)/du! should be regarded as the intrinsic tran- 
sition strength. In Sec. IV Bl we assume the strong cou- 
pling scheme |46j in order to transform calculated intrin- 
sic strength to the quantity in the laboratory frame. The 
response function is written as 

n(r,r';w)=no(r ) r';w) 

dr dr n (r, r ;uj)v(r , r )II(r , r';cj),(10) 
A 



n (r,r';a,) = ^{<Mr)G(-)(r,r'; 



to 



(r') 



i=i 



+#(r) GW(i 



■^)^(r')}- ( 



11) 



The single-particle Green's function in Eq. (|ll|l is defined 
by 



G^(r,r';E) = (r\(E-h[p }±ivr 1 \r n 



(12) 



Here, the superscript +(— ) indicates the outgoing (in- 
coming) boundary condition. In case that h[po] is rota- 
tionally invariant, the Green's function of Eq. (|12|) can 
be constructed by using the partial-wave expansion, 



G^{r,r';E)=2mJ2Y lm (v 



W[ui 



„(±)l 



(13) 

Here, W is the Wronskian and ui and wi are solutions of 
the radial Schrodinger equation for h[po] = — V 2 /2m + 
V(r): 



E- 



1 d 2 



2m dr 2 2mr 2 



1(1 + 1) 



V(r) Ri(r) = 0. (14) 



Ui is regular at origin and has an outgoing/incoming 
asymptotic form. In the ID continuum RPA [2l| . 



Eq. (|10[) is also expanded in partial waves. Then, the 
ID RPA response function (in the radial coordinate) is 
explicitly constructed by using Eqs. (UJJ, (JUJ, and lfl3")l . 

There are some difficulties to extend the theory to 
non-spherical systems. The first one is a purely numer- 
ical difficulty. Since the number of spatial grid points 
in the 3D space is much larger than that of the radial 
grid points, it is hard to explicitly construct the response 
function, n(r,r';w) and to perform spatial multi-fold in- 
tegration. We also need to calculate an inverse matrix 
to solve Eq. I|1U|) . The second difficulty lies in the com- 
plexity of boundary condition. Equations <|13[) and l|14l) 
cannot be used for cases of a deformed HF potential. 

We solve the first numerical problem by using an it- 
erative procedure for implicit calculation of the response 
and Green's function. For instance, in order to calculate 
the transition density, we recast Eq. © into an integral 
equation for Sp, 



dr |<5(r — r )— J rfr'n (r, r'; u)v(r', r )jSp(r ;w) 
dVn (r,r»F cxt (r'; W ). (15) 



This is equivalent to a linear algebraic equation in the 3D 
grid space and we use the iterative method to solve it. 
For the linear algebraic problem, A\x) = \b), the iterative 
methods require neither a full knowledge of the matrix A 
nor an inverse matrix A~ 1 , but do only results of operat- 
ing A on a certain vector \y). This means that we do not 
need to calculate an explicit form of Ho(r,r';ui). All we 
need to calculate is the action of no; i.e., no • v ■ dp and 
Ho • V ex t where the dot indicates the integral in Eq. i|15[l . 
This is an advantage of the iterative method over the di- 
rect method. In addition, the iterative method is known 
to be very efficient for a large sparse matrix. Then, the 
next task is to calculate action of no . According to Eq. 
(|ll|l . we have to operate G ± (E) on certain states \y). 
Now the problem is coupled to the second difficulty, that 
is the continuum boundary condition for deformed sys- 
tems. 

We start to divide the HF potential into a long-range 
spherical part and a short-range deformed one, h[po] = 
-V 2 /2m + Vo(r) + V"(r). In the present work, V (r) 
is taken as the Coulomb potential of a sphere of radius 
1.2 A 1 / 3 fm with a uniform change Ze. The single-particle 
Green's function for ho — — V 2 /2m + Vo(r) is constructed 

in the same way as Eq. which is denoted by (E) 
below. We have an identity for G, 

G^(r,r';E) = G^\r,r';E) 
+ J d\" G^\v,v" ■E)V(v")G {± \v'\v' ] E).{16) 

The boundary condition of G^ determines an asymp- 
totic behavior of G (± \ The action of G (± \ |x (± )) = 
G^'ly) for a given state \y), is obtained by solving a 
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linear algebraic equation 

{l-G^v}\ X ^)=G^\y). (17) 

Here, we use, again, the iterative method to solve this 
equation. 

In summary, to obtain the transition density, we solve 
Eq. Q15p. In order to do this, we need to calculate the 
operation of IIo , which then requires us to solve Eq. i|17|) 
with a proper boundary condition. The procedure re- 
sults in multiple-nested linear algebraic equations which 
are solved with iterative methods, such as the conju- 
gate gradient method. The detailed algorithm is given 
in Ref. Q. 



III. REAL-TIME TDHF+ABC 
A. Absorbing boundary condition (ABC) 

The TDHF equation can be efficiently solved in the 3D 
lattice space in real time [ill EtI |48| . The same technique 
has been applied to TDDFT of finite [43, |5(J and infinite 
electronic systems |5l|. In the real-time calculation, we 
propagate single-particle wave functions {6j . . . . t A us- 
ing the same technique as that in Ref. |47| . 

<f>i(r, t + At) = exp(-iAt • h[p(t + At/2)])<j}i{r, t), (18) 

where the exponential operator is expanded in a power 
series to (At) 4 . The time step in following applications is 
taken as At = 0.001 MeV -1 . There are many good rea- 
sons for solving the problem in real-time representation. 
First, the computation algorithm becomes very simple. 
In order to make the time evolution, only the opera- 
tion of the HF Hamiltonian on a certain single-particle 
state, h[p}\ip), is needed to be calculated, though the self- 
consistency between p(t) and h[p] brings slight complica- 
tion. Secondly, a single calculation of the time evolution 
provides information for a wide range of energy. Thus, 
for the strength function in a wide energy region, the 
real-time calculation is often more efficient than that in 
the energy domain. Last but not least, the TDHF wave 
packet in real space in real time gives an intuitive under- 
standing of dynamics. In pioneering works on heavy-ion 
collisions, the TDHF dynamics nicely demonstrated time 
evolution of nuclear inelastic scattering [Til ItH [H^ . 

The TDHF time evolution is relatively easy in the 
present computer power. The problem is how to impose 
the continuum boundary condition. The exact treatment 
of the continuum such as the Green's function method 
is very difficult (even impossible) in this case, because 
the energy of escaping particles cannot be determined 
uniquely. Thus, we attempt an approximate treatment. 
The usual approach to TDHF in real space is to as- 
sume the wave function to be zero at some distance R 
from the origin, which we call "Box boundary condition 
(BBC)" hereafter. Then, the time evolution must be 



completed before a significant portion of the wave reaches 
the boundary. Seeking higher accuracy, we must employ 
a larger R value, increasing the computation task. In- 
stead, in this paper, we employ the "Absorbing boundary 
condition (ABC)", which introduces a complex absorb- 
ing potential outside of the system. The method was 
first tested by Hamamoto and Mottelson for a schematic 
one-dimensional model of TDHF calculation 53] . Since 
then, however, its capability has not been fully examined 
in nuclear theory. On the other hand, in other fields of 
quantum physics, especially in atomic and molecular col- 
lision theories, the method has become one of standard 
methods for calculations of reactive scattering problems 
(see a recent review article [54| and references therein). 
We have demonstrated that the ABC is able to produce 
results identical to that of the exact continuum with the 
Green's function for TDDFT study of photoabsorption 
in molecules ^4|. In nuclear three-body reaction models 
[l5| . the method was also tested in detail for deuteron 
breakup reaction and provides an alternative method to 
the continuum-discretized-coupled channels (CDCC). A 
similar approach has been tested to calculate nuclear res- 
onance states in a simple model |55j . 

The success of the ABC approach is based on its sim- 
plicity. Actually, it requires only a minor modification 
of the real-time TDHF code, simply adding a complex 
potential, —ifj(r). We replace the HF Hamiltonian in 
Eq. O by 



h[p] — > h[p] - w?(r). 



(19) 



This prescription is equivalent to the use of Green's func- 
tion of Eq. i|12|) in which the infinitesimal imaginary part 
r/ is replaced by a finite and coordinate-dependent fj(r). 
The absorbing potential must be zero in a region where 
the ground-state density has a finite value, and it is fi- 
nite outside of the system. Of course, the addition of 
the complex potential violates the unitarity of time evo- 
lution. Thus, the norm of each single-particle state de- 
creases with time, which represents a physical process, 
the emission of particles. 

We adopt the same form of absorptive potential as pre- 
vious work s |l4l Il5l . This is a linear dependence on the 
coordinate I56ll57ll: 



fj(r) 



for < r < R, 

ii]o r 2T for i? < r < i? + Ar. 



(20) 



The size of the inner model space (r < R) is chosen so 
that the HF ground state converges within this space. 
The outer space of width Ar (R < r < R + Ar) is the 
absorbing region that should be large enough to prevent 
reflection of emitted outgoing waves. The condition of a 
good absorber for a particle with mass m and energy E 
is given by 



7 



E 1 / 2 1 
. 7= < hoi < —ArV^lE 3 / 2 . 
Arv°m tU 



(21) 



Here, we demand the reflection smaller than 0.1% and the 
transmission smaller than 3.3%. The similar condition 
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was given in Refs. [Til IBH l57| . Since the condition is 
energy dependent, we choose Ar and W as, 



12 fm, 770 = 10 MeV. 



(22) 



This satisfy the condition of Eq. J2U for 7 < E < 60 
MeV. 

For the linear response calculation, first, we solve 
the static HF problem with the imaginary-time method 
|58| to determine the occupied HF orbitals {0°}i=i,...,A- 
Then, an external perturbative field, T4 x t(r,t) = 
eF(r)5(t), is turned on instantaneously at t = 0. This 
results in an initial state of the TDHF calculation as 



l> i (r,t = 0+) = e- i * F(r V°(r), 



(23) 



where the constant e is arbitrary but should be small 
enough to validate the linear approximation of Eq. 
We calculate time evolution of the expectation value of 
F(r) (assumed to be real), 



A 

(¥(t)|JW)) = / dr^*(r,i)F(r)<Mr,i) 

i— 1 

= J drF(r)Sp(r,t). (24) 

Here, we assume that the ground-state expectation value 
of F(r) is zero at the last equation. Comparing its Fourier 
transform with Eq. Jj|J), we have 



dB(uj;F) 
doj 



1 



Im / ^(*(t)|F|*(t))e i 



(25) 



Note that (^(i) |F|^(t)) is fully determined by wave func- 
tions in the inner space (r < R), as far as the linear ap- 
proximation is valid. This can be easily understood using 
a relation, 5p = Y^j 4>i*84>i + h.c, and the condition that 
$ = at r > R. 



and the same form for y(y) and z{w) as well. This func- 
tion has an asymptotic values, x{u) ~ u at 11 C io and 
x{u) ~ ku at u 3> xq. All the derivatives and integrals 
in (x, y, z)-space are mapped to those in (u, v, u>)-space. 
For instance, 



dx 2 
dr 



duV 

dx 2 du \dx J du 2 
dx dy dz 



dudvdw- 



du dv dw 



(27) 
(28) 



The 3D (u, v, u>)-space is discretized in square mesh and 
finite-point formula in this space is applied to numerical 
differentiation. The curvilinear grid space employed in 
Sees . II Vl and IVl is shown in Fig.^ 
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x[fm] 
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FIG. 1: Adaptive grid in the (x, y)-plane for the coordinate 
transformation of Eq. 126H with xo = 8 fm, k — 5 and n — 2. 
The (u, v, w)-space is discretized in square mesh of 0.9 fm. 



B. Adaptive-coordinate 3D grid space 

There is a significant improvement of the computa- 
tional cost by reducing the number of grid points in the 
outer space. Although we take the outer model space 
roughly the same size as the inner one in the radial co- 
ordinate (R w Ar), its volume is considerably larger be- 
cause the volume element increases as r 2 dr. Therefore, 
the calculation of wave functions in the outer model space 
consumes most part of the computation time. However, 
since wave functions in the outer space is irrelevant for 
(^(i)lFl^I'(i)), the accurate description is not necessary 
there. Therefore, we use the adaptive curvilinear coor- 
dinate in the small-amplitude TDHF+ABC calculation 
|59| . The coordinate transformation we use in this paper 
is 



1 + (k — 1)u/(xq smh(u/xo)) n ' 



IV. ILLUSTRATIVE EXAMPLES: 
APPLICATION WITH THE BKN INTERACTION 

In this section, some illustrative examples are shown to 
demonstrate effects of continuum, validity of bound-state 
(L 2 ) approximation, and comparison between Green's 
function and ABC approach. We adopt the BKN inter- 
action used in Ref. [47]]. Note that, for this schematic in- 
teraction, the spin-isospin degeneracy is assumed all the 
time and the Coulomb potential acts on all orbitals with 
a charge e/2. The HF one-body Hamiltonian is given by 

h[p} = -^-y 2 + ^t p+^-t 3 p 2 + W Y + W c , (29) 
2m 4 16 

where the Yukawa potential, Wy, and Coulomb poten- 
tial, Wc, consist of their direct terms only. The param- 
eters are taken from Table. I of Ref. |47j . 
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FIG. 2: Mass octupole strengths as functions of excitation 
energy for ls O calculated with the BKN interaction. The top 
panels show results of the self-consistent continuum RPA. The 
middles and bottoms show results of the discretized RPA of 
R = 20 and 8 fm, respectively. The smoothing parameter F 
increases from left to right, 0.5, 2, and 10 MeV. The solid 
curves show the RPA strength function, while the dashed 
show unperturbed one. Since each orbital has a four-fold 
degeneracy (spin-isospin) , the E3 strengths are those shown 
multiplied by e 2 /4. 



A. Continuum in the spherical nucleus: 16 O 

1. L 2 -approximation of the continuum 

RPA calculations are often performed on the L 2 basis 
set, such as the harmonic oscillator basis. Bound ex- 
cited states are well described in those calculations, but 
how accurate is the ^-approximation for resonance and 
continuum states? In other words, what size of model 
space is necessary to describe an excited state with a fi- 
nite life time? In Ref . 01 > we gi ye a relation between the 
box size R and the energy resolution AE for continuum 
states; AE ~ hv/R, where v is the velocity of an escap- 
ing particle. If we consider a resonance with life-time of 
r, we should read AE ~ hv/(R + vt). For a long-lived 
state, vt 3> R, this is identical to the uncertainty princi- 
ple, AE ~ H/t. However, for a state of vt <C R, such as 
broad resonance and non-resonant continuum, the reso- 
lution is limited by the size of model space. 

Using the BKN interaction, it is easy to perform the 
self-consistent ID continuum RPA calculation for closed- 
shell spherical nuclei. We show, in Fig. [21 results of 
the ID continuum RPA and the RPA in a box radius 
R with BBC, which is referred to as "discretized RPA", 
for isoscalar (IS) octupole resonance in 16 O. We should 
note that a similar study on monopole resonance can 
be found in Rcf. For the continuum RPA cal- 



culation, the outgoing boundary condition is imposed 
at r — 8 fm (top panels). For simplicity, we use the 
free asymptotic form, tu, ~ e ±lkr with k 2 /2m — 
[E - V c (r) - 1(1 + l)/2mr 2 ] at r = 8 fm, instead of the 
exact Coulomb wave function. For the discretized RPA, 
the radius of model space is chosen as R = 20 fm (mid- 
dle panels) and R = 8 f m (bottom) . Since the discretized 
RPA produces only discrete peaks, we use a smoothing 
parameter T, adding an imaginary part, iT/2, to the real 
energy u>. In the continuum calculation, though we do 
not need to smear out the continuum strength, we use 
the same value of T to make the resolution as coarse as 
the discretized RPA. 

The continuum RPA calculation clearly shows the 
low-energy (LEOR) and high-energy octupole resonance 
(HEOR). The single-particle energy for p-shell is about 
— 16 MeV in this calculation. Thus, the LEOR is a bound 
peak whose width is entirely from a smoothing parameter 
T. As you see in Fig. |21 the bound LEOR peak depends 
neither on the boundary condition nor on the box size R. 
This justifies the use of the discretized RPA for bound ex- 
cited states. On the other hand, the structure of HEOR, 
which is embedded in the continuum, strongly depends 
on values of R. Note that the continuum RPA results 
with r = 0.5 MeV is almost identical to that with T = 0, 
which means that the width of HEOR is not artificial in 
contrast to the LEOR. The parameter T actually controls 
the energy resolution. For the discretized calculations 
with R — 20 fm, we need T > 8 MeV to produce roughly 
identical results to the continuum calculation. For those 
with R = 8 fm, we still see some discrepancy even with 
r = 10 MeV. In order to obtain sensible results in the 
discretized basis, we should average the strength func- 
tion with r inversely proportional to the box size. We 
find an empirical formula, T « Z(hv/R) for this calcula- 
tion. The continuum results with T = 0.5 MeV can be 
reproduced by the discretized calculation if we employ 
a model space of R > 200 fm. It is nearly impossible 
to treat this size of the 3D grid space with a present 
computer power. Therefore, it is certainly desirable to 
develop a method of treating the continuum boundary 
condition for deformed nuclei. The methods described in 
Sees. ILT1 and IIIII will serve this purpose. Next, we discuss 
applications of these methods to the BKN interaction. 



2. Small- amplitude TDHF+ABC vs. continuum RPA 

In this section, we examine accuracy and feasibility of 
methods in Sec. Inland in Sec. IIIII We compare results of 
two methods and show how accurate the TDHF+ABC 
can be. 

We use the same BKN model, O, and, again, cal- 
culate octupole states in 16 O. For the 3D continuum 
RPA calculation, the model space is the 3D coordinate 
space of R = 8 fm, discretized in square mesh of Ax = 
Ay = Az = 1 fm. For the real-time method with small- 
amplitude TDHF+ABC, we use the adaptive curvilinear 
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FIG. 3: The same as Fig. |21 but calculated with (a) the 
Green's function method in the 3D grid space and (b) the 
real-time small-amplitude TDHF+ABC in the adaptive 3D 
space. The smoothing parameter is F = 0.5 MeV. Calculated 
spurious dipole strength is shown by crosses for < id < 5 
MeV in units of fm 2 . 



coordinate of Fig. ^ with R = 8 fm and Ar = 12 fm, 
and ?7o = 10 MeV for the absorbing potential. The time 
evolution is calculated up to t — 30 MeV -1 , then we 
perform the Fourier transform of the timc-dependent oc- 
tupole moment with ^(r) = r 3 Y3o(f). Results of the 
Green's function method is shown in Fig. 0(a) and those 
of the real-time method in Fig. [3] (b) . These figures are 
almost identical to each other, but one may notice small 
difference. First, the strength at 20 < u> < 30 MeV 
is slightly higher for the real-time calculation. This is 
probably because the condition for the absorber, ijZIfr. 
breaks down for low-energy particles (E < 7 MeV). Sec- 
ondly, the peak position is higher for the real-time calcu- 
lation by about 0.3 MeV for LEOR and about 0.6 MeV 
for HEOR. This is due to the use of adaptive coordinate 
representation. The time evolution of octupole moment 
is shown in Fig. 01 with use of the square and the adap- 
tive coordinate. Discrepancy seen at t > 1 MeV -1 corre- 
sponds to 0.3 MeV difference in the LEOR energy. Com- 
paring results of these calculations with those of the ID 
continuum RPA in the radial coordinate, we see a very 
good agreement (see the top- left panel of Fig. |3J). This 
means that the nucleonic continuum states are properly 




FIG. 4: Time evolution of the octupole moment as a function 
of time for ie O. The calculation in the adaptive (square) mesh 
coordinate is shown by the solid (dashed) line. 



treated in both calculations of the 3D coordinate space; 
the Green's function method and the small-amplitude 
TDHF+ABC. There is a small peak in the 3D calculation 
at oj — 2.5 MeV. This is due to small admixture of the 
spurious translational mode. In Fig. 01(a), the strength 
calculated with V cx t = rYio is presented by crosses for 
< uj < 5 MeV in units of fm 2 . In the ID contin- 
uum RPA calculation with the partial-wave expansion, 
these octupole and dipole modes are separated. Thus, 
this mixing is not present in Fig. [21 However, in the 3D 
grid space, the translational and rotational invariance of 
the Hamiltonian is not exact. Adopting finer grid spac- 
ing diminishes the spurious peak height and moves its 
position toward zero energy. 

Figure 01 demonstrates an interesting feature in real 
time. The total energy is conserved within 300 keV up 
to t = 30 MeV -1 . With the BBC instead of the ABC, the 
energy conservation becomes even better. The absolute 
scale of its vertical axis does not have a significant mean- 
ing because it depends linearly on the arbitrary small 
parameter e. In the beginning, there is interference be- 
tween the LEOR and HEOR, however, for t > 1 MeV" 1 , 
only the LEOR mode survives. This feature clearly indi- 
cate stability of the bound collective excitation and de- 
cay of the collective mode in the continuum. The single- 
mode oscillation of the LEOR continues to the end of 
the time evolution (t = 30 MeV -1 ). The HEOR decays 
into the nucleon emission within time scale of t ~ 0.5 
h/MeV. Therefore, a part of the calculated energy width 
of HEOR, at least a few MeV, is associated with this 
nucleon escape width. 



B. Continuum in the deformed nucleus: Ne 

Now let us discuss a light deformed nucleus, 20 Ne. 
This illustrates usefulness and difficulties of the present 
approaches. Using the BKN interaction, 20 Ne has a 
superdeformed prolate shape with (3 ss 0.6. This nu- 
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cleus has a ground-state rotational band and the mea- 
sured B(E2\ 2 + — > + ) value is consistent with the 
deformation. Former calculations of the variation after 
parity projection have produced the Van-type octupole- 
deformed ground state with the BKN |6(| and with the 
Skyrme interaction |6lj . Since the system is deformed, 
the ID continuum RPA is no longer applicable. This is 
the first attempt of the 3D continuum RPA calculation 
for deformed nuclei. 

We use the same model space as the previous calcula- 
tion on 16 O. The IS monopole (r 2 ) and quadrupole field 
{r 2 Y2K (?)) are adopted as the external perturbations in 
Eq. QJ. Results of the 3D continuum RPA are shown 
in Fig. [5] The calculated single-particle energy of the 
last occupied orbital is —10.8 MeV. Thus, all the high- 
energy peaks in the figure are embedded in the contin- 
uum. The giant quadrupole resonance shows three peaks 
in order of K = 0, 1, and 2 in increasing energy (Fig. [5] 
(b)). Energy spacing between K = and 1 peaks is 
smaller than that between K = 1 and 2. This agrees 
with the simple scaling rule |62j |. The result also indi- 
cates no low-energy quadrupole vibration except for the 
zero-mode with K = 1. This is a characteristic feature 
in the superdeformation [63ll6^| . 

The monopole strength seems to consist of two com- 
ponents: a peak at 15 MeV and a broad hump in the 
energy region of E > 20. The peak position is lower than 
that of the unperturbed peak, by about 5 MeV. For the 
monopole strength in 16 O, calculated strength is shifted 
higher in energy with the BKN interaction 0|- There- 
fore, we consider this lowering in energy due to strong 
coupling to the quadrupole resonance. In fact, the peak 



lies at exactly the same energy as the K = quadrupole 
resonance (Fig. El (b)). We have reported a similar re- 
sult for the oblate nucleus, 12 C Although the BKN 
interaction may not be realistic for arguing real phenom- 
ena in 20 Ne, effects of such coupling in the continuum 
between different multipole resonances in deformed nu- 
clei would be an interesting sub ject in future. There are 
experimental data on this issue |66l 1671 168( . 

At the end of this section, we would like to men- 
tion a numerical problem of the real-time TDHF+ABC 
method. We have a difficulty to calculate a certain class 
of IS modes of excitation with the real-time method. This 
is associated with zero (Nambu-Goldstone) modes. For 
instance, calculating TDHF time evolution with the ex- 
ternal perturbative IS K = 0/K = 1 octupole field, the 
center of mass of the nucleus starts moving because of 
coupling to the translational motion. Of course, if we 
adopt a very small grid spacing, these modes are decou- 
pled, which is guaranteed by the self-consistent HF+RPA 
theory. In practice, we use a mesh of order of 1 fm in the 
3D Cartesian coordinates and a finer mesh size drastically 
increases a computational task. The problem is more se- 
rious in deformed cases than in spherical, because the 
angular momentum selection rule no longer works. In 
addition, the deformed nucleus has the rotational mode 
as another zero mode which is clearly seen in the K = 1 
mode in Fig. [3 In 16 0, we are able to perform the time 
evolution up to t > 30 MeV -1 , however, for the K = 
octupole mode in 20 Ne , t ~ 10 MeV -1 is a limit of time 
period in which the reliable calculation can be done. This 
is, of course, a matter of computational cost. If we do 
not use the adaptive curvilinear coordinate and adopt a 
larger space, we can carry out a stable calculation for a 
longer period. Because of this problem, we shall discuss 
applications with the Skyrme interaction to the isovector 
(IV) giant dipole resonances (GDR) in the next section, 
which is more stable and feasible. 



V. GDR STUDIED WITH SKYRME 
TDHF+ABC 

A. Effects of time-odd mean field in ls O 

The continuum RPA with the Skyrme energy func- 
tional is a standard method for describing collective ex- 
citations in closed-shell spherical nuclei. However, its 
fully-self-consistent calculations have not been achieved 
in practice, neglecting residual spin-orbit and Coulomb 
interactions. In addition, some of the time-odd densi- 
ties, which are known to be important for nuclear mo- 
ment of inertia and the local Galilean invariance |37| . 
are often neglected in the continuum RPA. In this sec- 
tion, we present an application of the small-amplitude 
TDHF+ABC method to the giant dipole resonance 
(GDR) in 16 O. Then, we compare the result with that 
of the former ID continuum RPA (which neglected the 
residual spin-orbit, Coulomb, and spin-spin parts), dis- 
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FIG. 6: Results of the Skyrme TDHF+ABC for GDR in 16 0. 
(a-1) Time evolution of the El moment as a function of time 
calculated with the Sill-full, (a-2) The same as (a-1) but 
with Sill-even, (b) Calculated photoabsorption cross section 
as a function of excitation energy. The Sill-full calculation 
(solid line) is compared to the Sill-even (dashed line). The 
smoothing parameter T = 0.5 MeV is used, (c) The same 
as (b) but neglecting the spin density s. Experimental pho- 
toneutron cross section (x3.5) is shown by the dotted line 
[gsj. (d) The same as (b) but neglecting the current density 
j. See text for details. 



cussing effects of the residual interactions. 

We adopt the Skyrme energy functional as same as 
Eqs. (A.2), (A.15), and (A.16) in Ref. [z3- The static 
HF+BCS code based on this functional is called EV8 
which assumes the parity and the z-signature symmetry. 
In the present work, we do not assume any symmetry, in 
order to allow a time-dependent state to be any Slater de- 



terminant during the time evolution. The energy density 
is written in terms of local densities as 



H(r) 



1 

2m 



r(r)+H even (r)+H odd (r), (30) 



with 



n cvcn (r) 
H odd (r) 



n cvcn [p n ,pAp,pr,pV ■ J], (31) 
2 ,s-Vxj]. (32) 



H odd [j 2 



Here, we follow the notation in Ref. 37]. According to 

Ref. [z3, terms of s • T - 7 2 , s • As, and (V • s) 2 are 
omitted. The energy functional, H cvcn + H odd , keeps the 
local gauge invariance [3(| |33 ■ It is customary in the 
static HF calculation to take account of the center-of- 
mass correction by multiplying the first term in Eq. i|3U|) 
by a factor (^4 — 1)/A. We use this correction both for 
static and dynamic calculations. 

In order to see effects of time-odd components, we 
adopt the SIII interaction which was used in the ID con- 
tinuum RPA calculation in Ref. [23j. We perform the 
TDHF calculation with the full functional of H oven + 
H° dd , and the one neglecting H odd . Hereafter, let us call 
the former functional "Sill-full" , and the latter "SIII- 
even" . The instantaneous external field is chosen as 



V at (T,t)=eM(El,n = 0)5{t) 



rY w (v)S(t), 



(33) 

where e' B1 ) indicates the El recoil charge, Ne/A for pro- 
tons and — Ze/A for neutrons, e is an arbitrary small 
number. Then, solving the TDHF equation 



d 



2m 



v + v cvcn (t) + v (t) > 4>i(t), 



(34) 

for i = 1, • • • , A. The time evolution is performed up to 
t = 30 MeV" 1 . The time-even mean field, V cvcn , has 
been well tested against a large number of experimental 
observations. In order to test the time-odd mean field, 
T/ odd , we neec j to investigate dynamical properties of nu- 
clei. 

Figure (a) shows time evolution of calculated El 
dipole moment, (^(t)\M(El)\^(t)) . In the Sill-even cal- 
culation, we see a beating pattern which results in two 
main peaks of the dashed line in Fig. [fj] (b) . This is in a 
good agreement with the result of the ID continuum RPA 
[23J. However, the inclusion of the time-odd mean field, 
which is necessary for the Galilean invariance, changes 
the strength distribution into a single peak (solid line). 
We decompose effects of time-odd densities into those 
of current density j and spin density s in Fig. [B] (c) and 
(d), respectively. The current density provides additional 
residual interaction to push the GDR to higher energy by 
0.5—1.3 MeV, while the spin density merges the two main 
peaks into one. This effect of time-odd density is not spe- 
cial to the SIII interaction. The same effect is observed 
with the SGII parameters of the Skyrme interaction. See 
Ref. [TJ f° r a brief report on the same calculation with 
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the SGII force. It is somewhat surprising that not only 
the current but also the spin density significantly modify 
the GDR structure. Photoneutron cross section data (^j 
are shown in the panel (c) by a dotted line. Their abso- 
lute values are multiplied by 3.5, since the data indicates 
less than 20 % of the TRK sum rule. The experimental 
shape of the GDR resembles that of the Sill-even cal- 
culation, but the two main peaks are calculated lower 
by about 3 MeV. Agreement on the main peak position 
is slightly improved in the Sill-full calculation, though 
the calculated peak is still lower than the experiment by 
2 - 2.5 MeV. 

If the interaction commutes with the El operator, the 
oscillator sum 

/>oo 

S(E1) = / EB(EI;0+ -> E(l~))dE, (35) 
J o 

1 POO 

= V / E\(E(ln)\M(El,»)\0)\ 2 dE(36) 



is identical to the TRK classical sum rule value 



S(E1) C 



9e 2 NZ 
87rm A 



(37) 



shell closure at N = 8. Since both Z — 4 and N — 10 are 
the magic numbers at the prolate superdeformed shape, 
we expect 14 Be to be deformed as large as 8 Be. A new 
mode of excitation of significant interest is the soft El 
mode near the neutron drip line |77ll78| . Coupling in the 
continuum between the soft El mode and the quadrupole 
deformation is an unsolved problem which can be ad- 
dressed by the present method to some extent. 

We use the Skyrme interaction of the SIII parame- 
ters including the time-odd components (SIII- full). The 
adopted model space is the adaptive grid in Fig. ^ with 
R = 10 fm and Ar = 12 fm. Usually, the static HF 
calculation is carried out with constraint on the center 
of mass at the origin. However, in this calculation, we 
do not impose any condition on the center-of-mass and 
on the direction of the principal axis. Although this re- 
sults in heavy computation for the imaginary-time step, 
it turns out that this is important for the stable time 
evolution of the TDHF state kicked off by the exter- 
nal perturbation. The external field is the same form 
as Eq. but includes rY 1±1 . The TDHF calculation 

with the perturbative El field provides the El intrinsic 
strength, dB(u,M(El,K))/dw th roug h Eq. |gSJ|. As- 
suming the strong coupling scheme |4g, the B(E1) tran- 
sition strength in the laboratory frame is given by 



du> 



For lb O, the classical sum rule gives S(El) class = 59.4 e 2 dB{u>; El) 
fm 2 MeV. Since the Skyrme interaction has momentum 
and isospin dependence, the classical sum rule is violated 
to a certain extent. We have S(E1) = 75.1 e 2 fm 2 MeV 
for Sill-full, and S{E1) = 67.1 c 2 fm 2 MeV for Sill-even. 
The enhancement of the TRK sum is 26 % for Sill-full 
and it reduces to 13 % if we neglect the time-odd mean 
field. This difference mainly comes from the spin density. 
If we integrate the strength in the energy region up to 30 
MeV, we have S(E1) w S(El) class for both Sill-full and 
Sill-even. 



J dE x B(E1;0 + -> E x (l-))6(u - E x ), 
= [ dE x \{E x \M(El,K)\0)\ 2 8(u - E x ), 
^ dB(uj;M(El,K)) 



(38) 



K=0.±1 



B. El resonances in even-even Be isotopes 

Finally, we apply the small-amplitude TDHF+ABC 
method to El resonances in beryllium isotopes. Beryl- 
lium nuclei have been extensively studied both theoret- 
ically and experimentally (see, e.g., Ref. [x3 an d refer- 
ences therein). 8 Be is well-known for the a-a clustering 
structure with an elongated prolate shape. Valence neu- 
trons added to 8 Be are expected to cause variety of struc- 
ture change in the ground and excited states O 0, EH . 
10 Be has two neutrons in addition to 8 Be. The a-a dis- 
tance in the ground state is considered to be slightly 
smaller than that in 8 Be. 12 Be is a semi-magic nucleus 
(V = 8), however its properties are different from spher- 
ical closed-shell nuclei. The measured spectroscopic fac- 
tors suggest that the last neutron pair is two-thirds in 
the sd configurations [7(1 . The neighboring odd nucleus, 
11 Be, is famous for the parity inversion and for the halo 
structure in the ground state. The existence of 14 Be at 
the drip line beyond N = 8 also indicates weakening of 



Here, the state |0) (\E X )) is the intrinsic ground (excited) 
state. 

The static HF calculation predicts all these nuclei to 
be deformed in prolate shape in the ground state. Calcu- 
lated quadrupole deformations are given in Table ID As 
we expected, 8 Be and 14 Be possess large deformation. 
12 Be has the smallest deformation, but its proton distri- 
bution has a moderate deformation. The static HF anal- 
ysis on Be isotopes with the SIII force have been already 
done in Ref. [7!j. The total binding energies are well 
reproduced. Calculated occupied single-particle energies 
are listed in Table ITT1 In the linear response approxima- 
tion of the TDHF, the neutron continuum plays its role at 
energies higher than the absolute value of single-particle 
energy of the last-occupied neutron. Since the proton 
orbitals become deeply bound in neutron-rich nuclei, the 
proton continuum is expected to be less important. How- 
ever, in these Be isotopes, the protons are important to 
produce prolate deformation of the mean field. 

Now let us discuss dynamical properties of these nu- 
clei. In Fig. [7| calculated time-dependent El moment 
is presented as functions of time. The time evolution is 
calculated up to t = 30 MeV" 1 . The be ginning third of 
the total period is shown in Fig. [7] Note that the am- 
plitude is magnified by a factor of 10 in the latter half 
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TABLE I: Calculated quadrupole deformation for even-even 
Be isotopes. The last two columns show deformation of neu- 
tron and proton density distribution, p n and p p , separately. 





P 


Pn 


Pp 


8 Be 


1.07 


1.06 


1.09 


10 Be 


0.39 


0.33 


0.49 


12 Be 


0.12 


0.07 


0.22 


14 Be 


0.74 


0.77 


0.66 



TABLE II: Calculated neutron (n) and proton (p) single- 
particle energies in units of MeV for Be isotopes. Each state 
has a two-fold degeneracy associated with the time-reversal 
symmetry. 



8 Be 


10 


Be 


12 


Be 


14 


Be 


n p 


n 


P 


n 


P 


n 


P 


-24.8 -22.9 


-26.3 


-31.6 


-25.5 


-36.7 


-24.9 


-39.1 


-13.3 -11.5 


-12.5 


-16.0 


-11.5 


-20.6 


-14.0 


-25.7 




-9.7 




-10.3 




-9.1 










-4.3 




-3.5 














-2.6 





of period. Performing the Fourier transform, calculated 
B(E1) transition strength is shown in Fig. [SJ In 8 Be, 
we observe large splitting of the GDR peak associated 
with the large quadrupole deformation (/3 ps 1). The 
magnitude of splitting is more than 10 MeV, that brings 
down the K = peak (the oscillation along the symme- 
try axis) to around 15 MeV in energy. In Fig. 0(a), the 
amplitude of the K = oscillation almost monotonically 
decays as time increases, while that of the K — 1 mode 
shows a beating pattern. This results in a splitting of 
the high-lying K — 1 peak. In the total B(E1) strength 
function, this is seen as a small peak in the middle of 
two main peaks. We also see that the K = oscilla- 
tion stays longer than the K = 1. This is because the 
peak position is near the particle decay threshold, thus, 
the allowed phase space is smaller for the K = peak 
(Table |nj. 

In 10 Be, we see a similar behavior to 8 Be. Because 
of the smaller deformation, the lower K — peak shifts 
to higher energy by about 5 MeV. Figure [5] (b) again 
indicates the K = 1 mode split into two peaks. Although 
the ground-state deformation in 10 Be is less than half of 
8 Be, the energy splitting between the lowest and highest 
peaks is still as large as 7 ~ 8 MeV. 

Next, let us discuss 12 Be. The calculated quadrupole 
deformation is the smallest among these even-even iso- 
topes. In contrast to 8,10 Be, we do not see a distin- 
guished double-peak structure. The GDR shows a peak 
at 21 MeV with a broader structure around 25 MeV. 
There exist low-energy El strength in the continuum be- 
low 10 MeV. The peak very near zero energy is due to 
small admixture of the translational mode. Because of 
this mixing, the response function suffers from spurious 
oscillatory behavior at energy below 2 MeV. Thus, we 
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FIG. 7: Calculated El moment as functions of time for (a) 
8 Be, (b) 10 Be, (c) 12 Be, and (d) 14 Be. The dashed (solid) line 
indicates the external field with K = (K — ±1). Scale of 
the vertical axis is arbitrary because it linearly depends on 
the small parameter e. It is magnified by a factor of five for 
the latter half of period, 5 < t < 10 MeV -1 . 



concentrate our focus on states at ui > 2 MeV. 

Though the B(E1) strength in low energy looks small 
compared to that in the main GDR, the integrated 
strength in the energy region of 2 < lj < 10 MeV amounts 
to B(E1) fts 0.14 e 2 fm 2 . The lowest sharp peak at 4.5 
MeV has B(E1;0+ -> 1") « 0.023 e 2 fm 2 . The next 
lowest peak at 5.6 MeV has B(E1) « 0.027 e 2 fm 2 . Both 
peaks have a dominant K = 1 character. The low-lying 
1~ state has been recently observed in 12 Be The 
observed excitation energy is E x = 2.68(3) MeV with 
B(E1; 0+ -> 1") = 0.051(13) e 2 fm 2 . Our result is higher 
in energy by a factor of two and the sum of B(E1) for 
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FIG. 8: Calculated values of dB(El;0 -> l~)/dw for 
8,io,i2,i4 Be The gjnootjjjjjg p arame ter F = 0.2 MeV is used. 
The thin dashed (solid) line is a contribution of K" = CP 
(1~) states and the thick solid line for the total strength. 



the lowest two peaks is comparable to the experiment. 
The two-neutron pairing model in Ref. |8l| predicted the 
1~ energy very well (2.7 MeV) but about five times over- 
estimated B(E\). The shell model calculation with ex- 
tended single-particle wave functions in Ref. jgj wen re ~ 
produced the lowest 1~ state [E x = 2.14 ~ 2.9 MeV with 
B(E1) = 0.063 ~ 0.072 e 2 fm 2 depending on the inter- 
action and model space). They also calculated B(E1) 
strength distribution in the GDR energy region without 
taking account of the continuum. Although their results 
strongly depend on the adopted interaction and model 
space, the calculated GDR energy is lower than ours. A 
striking difference from our result is that they have pre- 
dicted three main peaks with the WBP interaction. It 
is not clear at present whether this difference is due to 
the treatment of the continuum or to the ground-state 
correlation. 

Finally, let us move to the drip line, 14 Be. The doubly- 



magic closed-shell configuration (N = 10, Z = 4) at su- 
perdeformation leads to the large quadrupole deforma- 
tion of (3 = 0.74. The K — and K = 1 resonance peaks 
are at different positions whose centroids are at 15 MeV 
and 24 MeV. Figure0(d) indicates quick damping of the 
K = oscillation. The oscillating pattern almost disap- 
pears by t = 3 MeV -1 . This leads to the large width 
of the K = peak in Fig. [S] (d) . As a consequence of 
the large width, the double-peak structure in the total 
B(E1) strength function is not as clear as in 8,10 Be. It 
looks more like a single broad resonance at 20 MeV with 
the width of about 20 MeV. In Fig. 0(d), after the K = 
mode disappears, the K = 1 mode becomes dominant at 
t > 3 MeV -1 . This long-lived high-frequency K = 1 
mode results in sub-peaks embedded in the broad K = 1 
resonance (20 < u < 25 MeV). 

It is known that the weakly bound neutrons strongly 
couples to the continuum and produces the large dipole 
strength (83|. The Coulomb breakup of 11 Be is a typi- 
cal example [iij. This is often called "threshold effect" 
which has a peak at the threshold energy. Since the 
SIII interaction gives the last neutron binding of 2 — 3 
MeV, the threshold effect is weak. In the present cal- 
culation, we do not have significant threshold strength. 
On the other hand, another soft dipole peak is seen at 
5 MeV. This peak carries B(E1) m 0.26 e 2 fm 2 . A 
Coulomb dissociation experiment seems to suggest en- 
hanced strength at E x w 2 and 5 MeV 0. The shell- 
model calculation of Ref. 82] also indicates a similar peak 
(E x = 6.76 - 7.46 MeV with B(E1) = 0.097 - 0.146 e 2 
fm 2 ). Using the SGII interaction, this peak is at 7 MeV 
with B(El) w 0.14 e 2 fm 2 20], which well agrees with 
the shell-model result. On the GDR main peaks, our 
result looks rather different from the shell-model: the 
shell-model indicates a single main peak at 12 — 17 MeV, 
while we have a broad resonance whose centroid is around 
20 MeV. Since the shell model also indicates continuous 
El strength in the energy region of oj > 10 MeV, this 
difference may be simply due to lack of the continuum in 
Ref. [13. 

Calculated TRK sum rule values are listed in TableHTTI 
The enhancement is slightly smaller than the spherical 
16 case. Among the even-even Be isotopes, the enhance- 
ment is the biggest for 12 Be whose deformation is the 



smallest. According to the analysis on 16 in Sec._ 
about the half of this enhancement comes form effect of 
the time-odd spin density. The large deformation leads 
to a strong coupling to the if-quantum number, and this 
may restrict dynamics of spin degrees of freedom. The 
soft El strength, which is defined by the oscillator sum 
up to 15 MeV in the table, varies among these isotopes. 
The large value in 8 Be is due to the large ground-state 
deformation that brings the low-energy K — peak down 
close to 14 MeV. In 14 Be, it is the largest. This is due to 
combination of the deformation, the soft dipole peak at 
5 MeV, and the large width of the K — resonance at 
15 MeV. The deformation parameter 6 is estimated from 
the average energies of K — and K = 1 modes. We 
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TABLE III: Energy weighted sum rule values in units of 
e 2 fm 2 MeV. The second column shows values of the classical 
TRK formula. The small-amplitude TDHF+ABC calculation 
produces values in the third columns. The fourth column 
gives the soft El strength of the energy weighted sum, which 
is defined by uj < 15 MeV. The last column indicates the 
deformation parameter obtained by the splitting of the GDR 
peaks. 





5(-El) c i ass 


S(El) 


S(E1;E < 15MeV) 


S 


8 Be 


29.7 


34.0 


3.14 


0.43 


10 Be 


35.7 


42.8 


1.26 


0.21 


12 Be 


39.6 


48.2 


2.54 


0.05 


14 Be 


42.5 


52.2 


7.57 


0.35 



use Eq. (6-344) in Ref. . The S turns out to be much 
smaller than the deformation of the HF density distribu- 
tion, (3. The deformation derived from the GDR splitting 
is known to well agree with that from the E2 moment for 
actinide nuclei 46]. In light nuclei, the geometrical in- 
terpretation of the GDR frequencies may not be justified 
so well. 



VI. CONCLUSION 

We have developed the linear response theory in the 
continuum applicable to deformed systems. The ex- 
act treatment of the continuum is done by the iterative 
method for construction of the Green's function in the 3D 
Cartesian grid space (3D continuum RPA). The method 
is identical to the conventional ID continuum RPA in the 
spherical limit. At the same time, we have shown that 
the approximate but yet accurate treatment of the con- 
tinuum can be done by the absorbing boundary condition 
(ABC) approach. The small-amplitude TDHF+ABC 
method in the linear response regime is practically iden- 
tical to the 3D continuum RPA. Applications of these 
methods to the TDHF with the BKN interaction re- 
veals their usefulness and accuracy. The real-time TDHF 
method has a difficulty when we study excitation modes 
coupled to the zero modes. Since the method is fully self- 
consistent, the increase of model space (finer grid spac- 
ing) will solve the problem, though it requires heavier 
computation. 

Applications to systems with a realistic effective in- 
teraction have been performed with the small-amplitude 
Skyrme TDHF+ABC. The analysis on the GDR in 16 
suggests a significant contribution coming from the time- 
odd mean field which was often neglected in the ID 
continuum RPA. The peak structure in the continuum 
is affected by these residual interactions, especially by 
the spin density. Since the spin-dependent terms in the 



Skyrme energy functional, such as s 2 , s • As, and (V- s) 2 , 
are not linked to the time-even components by the lo- 
cal gauge invariance, the analysis may give a useful con- 
straint on these parts of the Skyrme functional. 

The coupling to the continuum becomes more impor- 
tant for weakly bound systems. We have studied the 
deformed continuum of the GDR in Be isotopes. The 
large deformation splitting of about 10 MeV is predicted 
for 8:14 Be. The K = main peak is significantly low- 
ered by the deformation to less than 15 MeV. The time 
evolution of the El moment indicates different damp- 
ing between 8 Be and neutron-rich Be isotopes, especially 
for the K = dipole mode. The soft dipole strength 
(E < 10 MeV) appears in 12 Be and 14 Be. Considering 
the fact that the SIII parameters were not determined by 
the isovector properties, we have a reasonable agreement 
with experiment on the low-energy 1~ state in 12 Be and 
14 Be. 

In this paper, we have studied only the IV GDR in 
neutron-rich nuclei, because of the numerical difficulty 
discussed above. The IS modes in neutron-rich deformed 
nuclei are also interesting to investigate. For instance, 
the octupole correlation in superdeformed 14 Be is ex- 
pected to be stronger than 8 Be. This is because the 
superdeformed magic numbers are classified into two cat- 
egory, and the N = 10 shell closure has a stronger oc- 
tupole correlation than the N = 4 jHJ • The small- 
amplitude TDHF+ABC may be a good method to see 
how the continuum affects this expectation, 

An important extension of the present approaches is 
the inclusion of pairing. Since the pairing plays an im- 
portant role in heavy nuclei, this is very desirable but a 
difficult task. In this respect, we should mention that the 
HFB-based continuum QRPA has been recently proposed 
by Matsuo [87| to take account of the continuum for both 
particle- hole (p-h) and particle-particle (p-p)/hole-hole 
(h-h) channels. The combination of the 3D continuum 
RPA and the continuum QRPA may produce a general 
theory to calculate excited states in the p-h, p-p, and h-h 
continuum for nuclei in the whole nuclear chart. 
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